############################################################################
# Replication Code for                                                     #
# "The Effects of Forced Versus Selective Exposure to Propaganda in China" #
# Political Science Research and Methods                                   #
# Author: Xiaoxiao SHEN                                                    #
# Date: July 6, 2025                                                       #
############################################################################



rm(list=ls())

# Install necessary packages if not already installed
if (!require("ggplot2")) install.packages("ggplot2")
if (!require("broom")) install.packages("broom")
if (!require("tidyverse")) install.packages("tidyverse")
if (!require("plyr")) install.packages("plyr")
if (!require("dplyr")) install.packages("dplyr")
if (!require("fuzzyjoin")) install.packages("fuzzyjoin")
if (!require("stringr")) install.packages("stringr")
if (!require("knitr")) install.packages("knitr")
if (!require("gridExtra")) install.packages("gridExtra")
if (!require("estimatr")) install.packages("estimatr")
if (!require("Matrix")) install.packages("Matrix")
if (!require("mediation")) install.packages("mediation")
if (!require("GGally")) install.packages("GGally")
if (!require("stargazer")) install.packages("stargazer")
if (!require("irr")) install.packages("irr")
if (!require("GK2011")) install.packages("GK2011")
if (!require("psych")) install.packages("psych")
if (!require("nFactors")) install.packages("nFactors")



### Step1: Load packages ###

library(ggplot2)
library(broom)
library(tidyverse)
library(plyr)
library(dplyr)
library(fuzzyjoin)
library(stringr)
library(knitr)
library(gridExtra)
library(estimatr)
library(Matrix)
library(mediation)
library(GGally)
library(stargazer)
library(irr)
library(GK2011)
library(psych)
library(nFactors)



# Set data directory --- Please replace "XXX" with the appropriate path
data_dir_survey1 <- "/XXX/Replication_Xiaoxiao Shen_PSRM/Output/Survey1"
data_dir_survey2 <- "/XXX/Replication_Xiaoxiao Shen_PSRM/Output/Survey2"
data_dir_survey3 <- "/XXX/Replication_Xiaoxiao Shen_PSRM/Output/Survey3"
data_dir_survey4 <- "/XXX/Replication_Xiaoxiao Shen_PSRM/Output/Survey4"
data_dir_survey5 <- "/XXX/Replication_Xiaoxiao Shen_PSRM/Output/Survey5"
data_dir_survey6 <- "/XXX/Replication_Xiaoxiao Shen_PSRM/Output/Survey6"
output_dir <- "/XXX/Replication_Xiaoxiao Shen_PSRM/Output/Combined"
if (!dir.exists(output_dir)) {
  dir.create(output_dir, recursive = TRUE)
}



### Step2: Load data ###

header_survey1 <- read.csv(file.path(data_dir_survey1, "8_data_for_combining.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_survey1 <- read.csv(file.path(data_dir_survey1, "8_data_for_combining.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_survey1) = header_survey1
rm(header_survey1) 
data_survey1 <- data_survey1[-1, -1]

header_survey2 <- read.csv(file.path(data_dir_survey2, "8_data_for_combining.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_survey2 <- read.csv(file.path(data_dir_survey2, "8_data_for_combining.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_survey2) = header_survey2
rm(header_survey2) 
data_survey2 <- data_survey2[-1, -1]

header_survey3 <- read.csv(file.path(data_dir_survey3, "8_data_for_combining.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_survey3 <- read.csv(file.path(data_dir_survey3, "8_data_for_combining.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_survey3) = header_survey3
rm(header_survey3) 
data_survey3 <- data_survey3[-1, -1]

header_survey4 <- read.csv(file.path(data_dir_survey4, "8_data_for_combining.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_survey4 <- read.csv(file.path(data_dir_survey4, "8_data_for_combining.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_survey4) = header_survey4
rm(header_survey4) 
data_survey4 <- data_survey4[-1, -1]

header_survey5 <- read.csv(file.path(data_dir_survey5, "8_data_for_combining.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_survey5 <- read.csv(file.path(data_dir_survey5, "8_data_for_combining.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_survey5) = header_survey5
rm(header_survey5) 
data_survey5 <- data_survey5[-1, -1]

header_survey6 <- read.csv(file.path(data_dir_survey6, "8_data_for_combining.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_survey6 <- read.csv(file.path(data_dir_survey6, "8_data_for_combining.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_survey6) = header_survey6
rm(header_survey6) 
data_survey6 <- data_survey6[-1, -1]

data = rbind.data.frame(data_survey1, data_survey2, data_survey3, data_survey4, data_survey5, data_survey6)



### Step3: Data Cleaning ###

# 3.1 Assign treatment
table(data$treatment) # force_polit:1220, force_science:1254, select_polit:1512, select_science:1179


# 3.2 Recode gender (1 = female)
data$gender = as.numeric(data$gender)
table(data$gender) # 0:2345, 1:2820


# 3.3 Recode age (in years)
# 1 for obs 30 or younger, 0 for obs older than 30
data$age = as.numeric(data$age)
table(data$age) # 0: 938, 1: 4227


# 3.4 Recode education
# 1 for college degree "7" or above, 0 for others
data$education = as.numeric(data$education)
table(data$education) # 0: 2496， 1: 2669


# 3.5 Recode urban hukou
# 0 for urban hukou, 0 for others
data$urban = as.numeric(data$urban)
table(data$urban) # 0: 2629, 1: 2536


# 3.6 Recode CCP
# 1 for CCP member, 0 for others
data$party = as.numeric(data$party)
table(data$party) # 0: 4479, 1:686


# 3.7 Recode income
# 1 for high income, 0 for others
data$income = as.numeric(data$income)
table(data$income) # 0:3629, 1: 1536


# 3.8 Recode knowledge
# Respondents' political knowledge is scored as follows: 
# 0 if both questions are answered incorrectly; 
# 0.5 if one is answered correctly; 
# and 1 if both are answered correctly.
data$knowledge = as.numeric(data$knowledge)
table(data$knowledge) # 0: 885， 0.5: 2115， 1: 2165


# 3.9 Recode marriage
# 0 = No (never married/divorced/other), 1 = Yes (married/widowed)
data$marriage = as.numeric(data$marriage)
table(data$marriage) # 0:3599, 1:1566


# 3.10 Recode interst
data$interest = as.numeric(data$interest)
table(data$interest) # 1: 92, 2: 772, 3: 1683, 4: 2162, 5: 456


# 3.11 Recode news_consumption
data$news_consumption = as.numeric(data$news_consumption)
table(data$news_consumption) # 1:973, 2: 711, 3: 1628, 4: 1046, 5:807


# 3.12 Recode occupation 
# Government agency is coded as 1 
# (1: government department, 
# 2: state-owned public institution [science, education, culture, health], 
# 3: state-owned enterprise); 
# non-government agency is coded as 0.
data$occupation = as.numeric(data$occupation)
table(data$occupation) # 0: 4133, 1:1032


# 3.13 Recode major
# social science为1，非social science为0
data$major = as.numeric(data$major)
table(data$major) # 0: 4415, 1:750


# 3.14 Recode ethnicity
data$ethnicity = as.numeric(data$ethnicity)
table(data$ethnicity) # 0: 314, 1: 4851


# 3.15 summary_check_RA
# Use package irr to calculate intercoder reliability，method: Cohen's Kappa
table(data$summary_check_RA1)
table(data$summary_check_RA2)
data$summary_check_RA1 <- as.numeric(data$summary_check_RA1)
data$summary_check_RA2 <- as.numeric(data$summary_check_RA2)

kappa_result <- kappa2(data.frame(data$summary_check_RA1, data$summary_check_RA2))


# 3.16 Summary statistics
selected_vars <- c("gender", "age", "education", "urban", "party", "income", "marriage", "occupation", "major", "ethnicity")

one_proportions <- sapply(data[selected_vars], function(x) mean(x == 1, na.rm = TRUE))

summary_stat <- data.frame(
  variable = names(one_proportions),
  proportion_of_1 = round(one_proportions, 3),
  stringsAsFactors = FALSE
)

summary_stat <- rbind(
  summary_stat,
  data.frame(
    variable = "N",
    proportion_of_1 = nrow(data)
  )
)

write.csv(summary_stat, file = file.path(output_dir, "3_summary_stat.csv"), row.names = FALSE)



### Step4: Balance Plot ###

data$forced_prop <- data$group=="1_forced_propaganda" 
data$forced_nonprop <- data$group=="2_forced_nonpropaganda" 
data$select <- data$group=="3_select" 

forced_prop <- summary(lm(forced_prop~gender+age+urban+education+income+party+marriage+knowledge+interest+news_consumption+ethnicity+major+occupation,data=data))
forced_nonprop <- summary(lm(forced_nonprop~gender+age+urban+education+income+party+marriage+knowledge+interest+news_consumption+ethnicity+major+occupation,data=data))
select <- summary(lm(select~gender+age+urban+education+income+party+marriage+knowledge+interest+news_consumption+ethnicity+major+occupation,data=data))

forced_prop.1 <- as.data.frame(forced_prop$coefficients)
forced_prop.1$Group <- 'forced_prop'

forced_nonprop.1 <- as.data.frame(forced_nonprop$coefficients)
forced_nonprop.1$Group <- 'forced_nonprop'

select.1 <- as.data.frame(select$coefficients)
select.1$Group <- 'select'

all_4.rap = rbind.data.frame(forced_prop.1, forced_nonprop.1, select.1)
all_4.rap$Estimate = round(all_4.rap$Estimate, 3)
all_4.rap$`Std. Error` = round(all_4.rap$`Std. Error`, 3)
all_4.rap$`t value` = round(all_4.rap$`t value`, 3)
all_4.rap$`Pr(>|t|)` = round(all_4.rap$`Pr(>|t|)`, 3)
write.csv(all_4.rap, file = file.path(output_dir, "4_balance_check.csv"), row.names = TRUE)



### Step5: Manipulation check ###

table(data$class) # forced: 2474, select: 2691

# effort
data$effort = as.numeric(data$effort)
data$class = as.factor(data$class)

summary(aov(data$effort~class,data=data))
TukeyHSD(aov(data$effort~class,data=data))

effort.1=TukeyHSD(aov(data$effort~class,data=data))
effort.1=as.data.frame(effort.1$class)


# trivialization 
data$trivialization = as.numeric(data$trivialization)
data$class = as.factor(data$class)

summary(aov(data$trivialization~class,data=data)) 
TukeyHSD(aov(data$trivialization~class,data=data))

trivialization.1=TukeyHSD(aov(data$trivialization~class,data=data))
trivialization.1=as.data.frame(trivialization.1$class)

manipulation_check = rbind.data.frame(effort.1, trivialization.1)
manipulation_check$SE = (manipulation_check$upr - manipulation_check$lwr)/1.96/2
manipulation_check = round(manipulation_check, 3)
manipulation_check$Variable = c("Effort", "Trivialization")

write.csv(manipulation_check, file = file.path(output_dir, "5_manipulation_check.csv"), row.names = TRUE)



### Step6: Marginal effect ###

# 6.1 reverse function
reverse=function(x,old,new,store=-1){ 
  # the function reverses the original score to a new one
  # "store" is a temp store which keeps the dt, make sure "can" is not equal to any one of "old" and "new"
  # make sure old and new are the same length
  # "x" is a vector
  numb=length(old)
  for(i in 1:numb){
    x[which(x==old[i])] = store
    x[which(x==new[i])] = old[i]
    x[which(x==store)] = new[i]
  }
  return(x) # return a vector
}


# 6.2 rescale covariates
range_0_1=function(x,range=NULL){ 
  # if do not give range, make sure both best and worst answer has sample
  if(is.null(range)){
    min=min(x,na.rm = TRUE)
    max=max(x,na.rm = TRUE)
  }else{
    min=min(range)
    max=max(range)
  }
  if(min>0){
    print(min)
    warning("minimum value larger than 0, make sure it is correct")
  }
  result=(x-min)/(max-min)
  return(result)
}


# 6.3 select_polit=1, select_science=0
data$select_polit=NA # 830 obs. of 131 variables
data$select_polit[which(data$treatment=="select_polit")]=1
data$select_polit[which(data$treatment=="select_science")]=0 # no not need to choose extra variable, just use "select" is enough
table(data$select_polit) # 0: 1179, 1: 1512


# 6.4 prepare: political_ideology / political_interest / political_knowledge / news_consumption

# 6.4.1 political_ideology
data$ideology=as.numeric(data$ideology) 
table(data$ideology)


# 6.4.2 political_interest
data$interest = as.numeric(data$interest)
table(data$interest) # 0:92, 0.25:772, 0.5:1683, 0.75:2162, 1:456


# 6.4.3 political_knowledge
data$knowledge = as.numeric(data$knowledge)
table(data$knowledge) # 0:885, 0.5:2115, 1:2165


# 6.4.4 news_consumption
data$news_consumption <- as.numeric(data$news_consumption)
table(data$news_consumption) # 0:973, 0.25:711, 0.5:1628, 0.75:1046, 1:807


# 6.5 logistic regression

# 6.5.1 binomial links to logit regression
marginal.reg=glm(select_polit~gender+age+urban+education+income+party+marriage+knowledge+interest+news_consumption+ethnicity+major+occupation,data=data,family="binomial")


# 6.5.2 summary 
marginal_reg_summary = as.data.frame(summary(marginal.reg)$coefficients)
marginal_reg_summary = round(marginal_reg_summary, 3)
write.csv(marginal_reg_summary, file = file.path(output_dir, "8_marginal_effects.csv"), row.names = TRUE)



### Step9: PCA Analysis ###

# 9.1 create new dataframe
data$central_trust = as.numeric(data$central_trust)
data$central_satisfaction = as.numeric(data$central_satisfaction)
data$local_trust = as.numeric(data$local_trust)
data$local_satisfaction = as.numeric(data$local_satisfaction)
data$evaluation_now = as.numeric(data$evaluation_now)
data$evaluation_future = as.numeric(data$evaluation_future)
data$confidence_1 = as.numeric(data$confidence_1)
data$confidence_2 = as.numeric(data$confidence_2)
data$efficacy_1 = as.numeric(data$efficacy_1)
data$efficacy_2 = as.numeric(data$efficacy_2)
data$pride_1 = as.numeric(data$pride_1)
data$pride_2 = as.numeric(data$pride_2)
data$rand = as.numeric(data$rand)
data$tr = as.numeric(data$tr)
data$treatment = as.numeric(data$treatment)
data$emotion_1 = as.numeric(data$emotion_1)
data$emotion_2 = as.numeric(data$emotion_2)
data$emotion_3 = as.numeric(data$emotion_3)
data$emotion_4 = as.numeric(data$emotion_4)
data$emotion_5 = as.numeric(data$emotion_5)
data$emotion_6 = as.numeric(data$emotion_6)
data$effort = as.numeric(data$effort)
data$trivialization = as.numeric(data$trivialization)

dt <- as.data.frame(matrix(nrow=dim(data)[1]*2,ncol=0)) 
dt$central <- append(data$central_trust, data$central_satisfaction)
dt$local <- append(data$local_trust, data$local_satisfaction)
dt$evaluation <- append(data$evaluation_now, data$evaluation_future)
dt$confidence <- append(data$confidence_1, data$confidence_2)
dt$efficacy <- append(data$efficacy_1, data$efficacy_2)
dt$pride <- append(data$pride_1, data$pride_2)
dt$rand <- append(data$rand, data$rand)
dt$tr <- append(data$tr, data$tr)
dt$treatment <- append(data$treatment, data$treatment)
dt$emotion_1 <- append(data$emotion_1, data$emotion_1)
dt$emotion_2 <- append(data$emotion_2, data$emotion_2)
dt$emotion_3 <- append(data$emotion_3, data$emotion_3)
dt$emotion_4 <- append(data$emotion_4, data$emotion_4)
dt$emotion_5 <- append(data$emotion_5, data$emotion_5)
dt$emotion_6 <- append(data$emotion_6, data$emotion_6)
dt$effort <- append(data$effort, data$effort)
dt$trivialization <- append(data$trivialization, data$trivialization)


# 9.2 data pca analysis
result_dt=principal(dt[,c(1,2,4,5)], nfactors = 1, residuals = FALSE,rotate="varimax")
result_dt
dt$pca_score=result_dt$score[,1]
dt$pca_score

ev <- eigen(cor(dt[,c(1,2,4,5)])) 
nS <- nScree(x=ev$values)
plotnScree(nS)


# Proportion of Variance Explained
sink(file.path(output_dir, "9.0.1 PCA_Proportion of Variance.txt"))
pca_proportion <- princomp(cor(dt[, c(1, 2, 4, 5)]))
summary(pca_proportion)
sink()


# Loadings
sink(file.path(output_dir, "9.0.2 PCA_Loadings.txt"))
pca_proportion <- princomp(cor(dt[,c(1,2,4,5)]))
summary(pca_proportion)[2]
sink()


# Principal Component Analysis on Pro-regime Attitudes
gk_9.estimate=list()
gk_9.estimate[[1]] = estimate(rand = dt$rand, tr = dt$tr, y = dt$pca_score)
list_.name=c("PCA Score")
names(gk_9.estimate)=list_.name
gk_9.estimate_pca_score = as.data.frame(gk_9.estimate$`PCA Score`)
gk_9.estimate_pca_score['name'] = rep("PCA Score", times=4)
gk_9.estimate_pca_score
write.csv(gk_9.estimate_pca_score, file = file.path(output_dir, "9_pca_on_pro_regime_attitudes.csv"), row.names = TRUE)



### Step10: Potential Mechanisms ###

# 10.1 rand
table(data$rand) # 0: 2691, 1: 2474


# 10.2 tr
table(data$tr) # 0: 2433, 1: 2732


# 10.3 gk.estimate
set.seed(12345)

data$effort = as.numeric(data$effort)
data$trivialization = as.numeric(data$trivialization)
data$`emotion_1` = as.numeric(data$`emotion_1`)
data$`emotion_2` = as.numeric(data$`emotion_2`)
data$`emotion_3` = as.numeric(data$`emotion_3`)
data$`emotion_4` = as.numeric(data$`emotion_4`)
data$`emotion_5` = as.numeric(data$`emotion_5`)
data$`emotion_6` = as.numeric(data$`emotion_6`)

gk.estimate=list()
gk.estimate[[1]]=estimate(rand = data$rand, tr = data$tr, y = data$effort)
gk.estimate[[2]]=estimate(rand = data$rand, tr = data$tr, y = data$trivialization)
gk.estimate[[3]]=estimate(rand = data$rand, tr = data$tr, y = data$`emotion_1`)
gk.estimate[[4]]=estimate(rand = data$rand, tr = data$tr, y = data$`emotion_2`)
gk.estimate[[5]]=estimate(rand = data$rand, tr = data$tr, y = data$`emotion_3`)
gk.estimate[[6]]=estimate(rand = data$rand, tr = data$tr, y = data$`emotion_4`)
gk.estimate[[7]]=estimate(rand = data$rand, tr = data$tr, y = data$`emotion_5`)
gk.estimate[[8]]=estimate(rand = data$rand, tr = data$tr, y = data$`emotion_6`)

list.name=c("Effort","Trivialization",
            "Anxious","Proud",
            "Angry","Hope",
            "Worry","Excited")

names(gk.estimate)=list.name # remember to give every estimate a name

gk.estimate_Effort = as.data.frame(gk.estimate$Effort)
gk.estimate_Trivialization = as.data.frame(gk.estimate$Trivialization)
gk.estimate_Anxious = as.data.frame(gk.estimate$Anxious)
gk.estimate_Proud = as.data.frame(gk.estimate$Proud)
gk.estimate_Angry = as.data.frame(gk.estimate$Angry)
gk.estimate_Hope = as.data.frame(gk.estimate$Hope)
gk.estimate_Worry = as.data.frame(gk.estimate$Worry)
gk.estimate_Excited = as.data.frame(gk.estimate$Excited)

gk.estimate_Effort['name'] = rep("Effort", times=4)
gk.estimate_Trivialization['name'] = rep("Trivialization", times=4)
gk.estimate_Anxious['name'] = rep("Anxious", times=4)
gk.estimate_Proud['name'] = rep("Proud", times=4)
gk.estimate_Angry['name'] = rep("Angry", times=4)
gk.estimate_Hope['name'] = rep("Hope", times=4)
gk.estimate_Worry['name'] = rep("Worry", times=4)
gk.estimate_Excited['name'] = rep("Excited", times=4)

gk.estimate.rap=rbind.data.frame(gk.estimate_Effort, gk.estimate_Trivialization,
                                 gk.estimate_Anxious, gk.estimate_Proud,
                                 gk.estimate_Angry, gk.estimate_Hope,
                                 gk.estimate_Worry, gk.estimate_Excited)

gk.estimate.rap$Estimate = round(gk.estimate.rap$Estimate, 3)
gk.estimate.rap$SE = round(gk.estimate.rap$SE, 3)
gk.estimate.rap$t = round(gk.estimate.rap$t, 3)
gk.estimate.rap$p = round(gk.estimate.rap$p, 3)

# 10.3 gk.estimate.rap
write.csv(gk.estimate.rap, file = file.path(output_dir, "10_potential_mechanisms.csv"), row.names = TRUE)



### 11 High/Low -> ATE/selector/non-selector (Y is PCA) ###

# 11.1 political ideology
std=function(x){
  return((x-mean(x))/sd(x))
}

data$ideology = as.numeric(data$ideology)
table(data$ideology)

data$high_ideology=ifelse(std(data$ideology)>=0,1,0)
table(data$high_ideology) # 0: 2429, 1: 2736

high.ideology <- data[data$high_ideology == 1,]
low.ideology <- data[data$high_ideology == 0,]


# 11.1.1 high.ideology
dt_high.ideology <- as.data.frame(matrix(nrow=nrow(high.ideology)*2,ncol=0)) 
dt_high.ideology$central <- append(high.ideology$central_trust, high.ideology$central_satisfaction)
dt_high.ideology$local <- append(high.ideology$local_trust, high.ideology$local_satisfaction)
dt_high.ideology$evaluation <- append(high.ideology$evaluation_now, high.ideology$evaluation_future)
dt_high.ideology$confidence <- append(high.ideology$confidence_1, high.ideology$confidence_2)
dt_high.ideology$efficacy <- append(high.ideology$efficacy_1, high.ideology$efficacy_2)
dt_high.ideology$pride <- append(high.ideology$pride_1, high.ideology$pride_2)
dt_high.ideology$rand <- append(high.ideology$rand, high.ideology$rand)
dt_high.ideology$tr <- append(high.ideology$tr, high.ideology$tr)

result_dt_high.ideology=principal(dt_high.ideology[,c(1,2,4,5)], nfactors = 1, residuals = FALSE,rotate="varimax")
result_dt_high.ideology
dt_high.ideology$pca_score=result_dt_high.ideology$score[,1]
dt_high.ideology$pca_score

gk_1.1.1.estimate=list()
gk_1.1.1.estimate[[1]] = estimate(rand = dt_high.ideology$rand, tr = dt_high.ideology$tr, y = dt_high.ideology$pca_score)

list_.name=c("PCA Score")
names(gk_1.1.1.estimate)=list_.name
gk_11.1.1.estimate_pca_score = as.data.frame(gk_1.1.1.estimate$`PCA Score`)
gk_11.1.1.estimate_pca_score['name'] = rep("PCA Score", times=4)
gk_11.1.1.estimate_pca_score$Estimate = round(gk_11.1.1.estimate_pca_score$Estimate, 3)
gk_11.1.1.estimate_pca_score$SE = round(gk_11.1.1.estimate_pca_score$SE, 3)
gk_11.1.1.estimate_pca_score$t = round(gk_11.1.1.estimate_pca_score$t, 3)
gk_11.1.1.estimate_pca_score$p = round(gk_11.1.1.estimate_pca_score$p, 3)
write.csv(gk_11.1.1.estimate_pca_score, file = file.path(output_dir, "11_1_pca_on_pro_regime_attitudes_high_ideology.csv"), row.names = TRUE)



# 11.1.2 low.ideology
dt_low.ideology <- as.data.frame(matrix(nrow=nrow(low.ideology)*2,ncol=0)) 
dt_low.ideology$central <- append(low.ideology$central_trust, low.ideology$central_satisfaction)
dt_low.ideology$local <- append(low.ideology$local_trust, low.ideology$local_satisfaction)
dt_low.ideology$evaluation <- append(low.ideology$evaluation_now, low.ideology$evaluation_future)
dt_low.ideology$confidence <- append(low.ideology$confidence_1, low.ideology$confidence_2)
dt_low.ideology$efficacy <- append(low.ideology$efficacy_1, low.ideology$efficacy_2)
dt_low.ideology$pride <- append(low.ideology$pride_1, low.ideology$pride_2)
dt_low.ideology$rand <- append(low.ideology$rand, low.ideology$rand)
dt_low.ideology$tr <- append(low.ideology$tr, low.ideology$tr)

result_dt_low.ideology=principal(dt_low.ideology[,c(1,2,4,5)], nfactors = 1, residuals = FALSE,rotate="varimax")
result_dt_low.ideology
dt_low.ideology$pca_score=result_dt_low.ideology$score[,1]
dt_low.ideology$pca_score

gk_11.1.2.estimate=list()
gk_11.1.2.estimate[[1]] = estimate(rand = dt_low.ideology$rand, tr = dt_low.ideology$tr, y = dt_low.ideology$pca_score)

list_.name=c("PCA Score")
names(gk_11.1.2.estimate)=list_.name
gk_11.1.2.estimate_pca_score = as.data.frame(gk_11.1.2.estimate$`PCA Score`)
gk_11.1.2.estimate_pca_score['name'] = rep("PCA Score", times=4)
gk_11.1.2.estimate_pca_score$Estimate = round(gk_11.1.2.estimate_pca_score$Estimate, 3)
gk_11.1.2.estimate_pca_score$SE = round(gk_11.1.2.estimate_pca_score$SE, 3)
gk_11.1.2.estimate_pca_score$t = round(gk_11.1.2.estimate_pca_score$t, 3)
gk_11.1.2.estimate_pca_score$p = round(gk_11.1.2.estimate_pca_score$p, 3)
write.csv(gk_11.1.2.estimate_pca_score, file = file.path(output_dir, "11_1_pca_on_pro_regime_attitudes_low_ideology.csv"), row.names = TRUE)



# 11.2 political interest
data$interest = as.numeric(data$interest)
table(data$interest) # 0:92, 0.25:772, 0.5:1683, 0.75:2162, 1:456

data$high_interest=ifelse(std(data$interest)>=0,1,0)
table(data$high_interest) # 0:2547, 1:2618

high.interest <- data[data$high_interest == 1,]
low.interest <- data[data$high_interest == 0,]


# 11.2.1 high.interest
dt_high.interest <- as.data.frame(matrix(nrow=nrow(high.interest)*2,ncol=0))
dt_high.interest$central <- append(high.interest$central_trust, high.interest$central_satisfaction)
dt_high.interest$local <- append(high.interest$local_trust, high.interest$local_satisfaction)
dt_high.interest$evaluation <- append(high.interest$evaluation_now, high.interest$evaluation_future)
dt_high.interest$confidence <- append(high.interest$confidence_1, high.interest$confidence_2)
dt_high.interest$efficacy <- append(high.interest$efficacy_1, high.interest$efficacy_2)
dt_high.interest$pride <- append(high.interest$pride_1, high.interest$pride_2)
dt_high.interest$rand <- append(high.interest$rand, high.interest$rand)
dt_high.interest$tr <- append(high.interest$tr, high.interest$tr)

result_dt_high.interest=principal(dt_high.interest[,c(1,2,4,5)], nfactors = 1, residuals = FALSE,rotate="varimax")
result_dt_high.interest
dt_high.interest$pca_score=result_dt_high.interest$score[,1]
dt_high.interest$pca_score

gk_11.2.1.estimate=list()
gk_11.2.1.estimate[[1]] = estimate(rand = dt_high.interest$rand, tr = dt_high.interest$tr, y = dt_high.interest$pca_score)

list_.name=c("PCA Score")
names(gk_11.2.1.estimate)=list_.name
gk_11.2.1.estimate_pca_score = as.data.frame(gk_11.2.1.estimate$`PCA Score`)
gk_11.2.1.estimate_pca_score['name'] = rep("PCA Score", times=4)

gk_11.2.1.estimate_pca_score$Estimate = round(gk_11.2.1.estimate_pca_score$Estimate, 3)
gk_11.2.1.estimate_pca_score$SE = round(gk_11.2.1.estimate_pca_score$SE, 3)
gk_11.2.1.estimate_pca_score$t = round(gk_11.2.1.estimate_pca_score$t, 3)
gk_11.2.1.estimate_pca_score$p = round(gk_11.2.1.estimate_pca_score$p, 3)

write.csv(gk_11.2.1.estimate_pca_score, file = file.path(output_dir, "11_2_pca_on_pro_regime_attitudes_high_interest.csv"), row.names = TRUE)



# 11.2.2 low.interest
dt_low.interest <- as.data.frame(matrix(nrow=nrow(low.interest)*2,ncol=0))
dt_low.interest$central <- append(low.interest$central_trust, low.interest$central_satisfaction)
dt_low.interest$local <- append(low.interest$local_trust, low.interest$local_satisfaction)
dt_low.interest$evaluation <- append(low.interest$evaluation_now, low.interest$evaluation_future)
dt_low.interest$confidence <- append(low.interest$confidence_1, low.interest$confidence_2)
dt_low.interest$efficacy <- append(low.interest$efficacy_1, low.interest$efficacy_2)
dt_low.interest$pride <- append(low.interest$pride_1, low.interest$pride_2)
dt_low.interest$rand <- append(low.interest$rand, low.interest$rand)
dt_low.interest$tr <- append(low.interest$tr, low.interest$tr)

result_dt_low.interest=principal(dt_low.interest[,c(1,2,4,5)], nfactors = 1, residuals = FALSE,rotate="varimax")
result_dt_low.interest
dt_low.interest$pca_score=result_dt_low.interest$score[,1]
dt_low.interest$pca_score

gk_11.2.2.estimate=list()
gk_11.2.2.estimate[[1]] = estimate(rand = dt_low.interest$rand, tr = dt_low.interest$tr, y = dt_low.interest$pca_score)

list_.name=c("PCA Score")

names(gk_11.2.2.estimate)=list_.name

gk_11.2.2.estimate_pca_score = as.data.frame(gk_11.2.2.estimate$`PCA Score`)

gk_11.2.2.estimate_pca_score['name'] = rep("PCA Score", times=4)

gk_11.2.2.estimate_pca_score$Estimate = round(gk_11.2.2.estimate_pca_score$Estimate, 3)
gk_11.2.2.estimate_pca_score$SE = round(gk_11.2.2.estimate_pca_score$SE, 3)
gk_11.2.2.estimate_pca_score$t = round(gk_11.2.2.estimate_pca_score$t, 3)
gk_11.2.2.estimate_pca_score$p = round(gk_11.2.2.estimate_pca_score$p, 3)

write.csv(gk_11.2.2.estimate_pca_score, file = file.path(output_dir, "11_2_pca_on_pro_regime_attitudes_low_interest.csv"), row.names = TRUE)



# 11.2.3 high.interest - potential mechanisms
set.seed(12345)

high.interest$effort = as.numeric(high.interest$effort)
high.interest$trivialization = as.numeric(high.interest$trivialization)
high.interest$`emotion_1` = as.numeric(high.interest$`emotion_1`)
high.interest$`emotion_2` = as.numeric(high.interest$`emotion_2`)
high.interest$`emotion_3` = as.numeric(high.interest$`emotion_3`)
high.interest$`emotion_4` = as.numeric(high.interest$`emotion_4`)
high.interest$`emotion_5` = as.numeric(high.interest$`emotion_5`)
high.interest$`emotion_6` = as.numeric(high.interest$`emotion_6`)

gk.estimate=list()
gk.estimate[[1]]=estimate(rand = high.interest$rand, tr = high.interest$tr, y = high.interest$effort)
gk.estimate[[2]]=estimate(rand = high.interest$rand, tr = high.interest$tr, y = high.interest$trivialization)
gk.estimate[[3]]=estimate(rand = high.interest$rand, tr = high.interest$tr, y = high.interest$`emotion_1`)
gk.estimate[[4]]=estimate(rand = high.interest$rand, tr = high.interest$tr, y = high.interest$`emotion_2`)
gk.estimate[[5]]=estimate(rand = high.interest$rand, tr = high.interest$tr, y = high.interest$`emotion_3`)
gk.estimate[[6]]=estimate(rand = high.interest$rand, tr = high.interest$tr, y = high.interest$`emotion_4`)
gk.estimate[[7]]=estimate(rand = high.interest$rand, tr = high.interest$tr, y = high.interest$`emotion_5`)
gk.estimate[[8]]=estimate(rand = high.interest$rand, tr = high.interest$tr, y = high.interest$`emotion_6`)

list.name=c("Effort","Trivialization",
            "Anxious","Proud",
            "Angry","Hope",
            "Worry","Excited")

names(gk.estimate)=list.name # remember to give every estimate a name

gk.estimate_Effort = as.data.frame(gk.estimate$Effort)
gk.estimate_Trivialization = as.data.frame(gk.estimate$Trivialization)
gk.estimate_Anxious = as.data.frame(gk.estimate$Anxious)
gk.estimate_Proud = as.data.frame(gk.estimate$Proud)
gk.estimate_Angry = as.data.frame(gk.estimate$Angry)
gk.estimate_Hope = as.data.frame(gk.estimate$Hope)
gk.estimate_Worry = as.data.frame(gk.estimate$Worry)
gk.estimate_Excited = as.data.frame(gk.estimate$Excited)

gk.estimate_Effort['name'] = rep("Effort", times=4)
gk.estimate_Trivialization['name'] = rep("Trivialization", times=4)
gk.estimate_Anxious['name'] = rep("Anxious", times=4)
gk.estimate_Proud['name'] = rep("Proud", times=4)
gk.estimate_Angry['name'] = rep("Angry", times=4)
gk.estimate_Hope['name'] = rep("Hope", times=4)
gk.estimate_Worry['name'] = rep("Worry", times=4)
gk.estimate_Excited['name'] = rep("Excited", times=4)

gk.estimate.rap=rbind.data.frame(gk.estimate_Effort, gk.estimate_Trivialization,
                                 gk.estimate_Anxious, gk.estimate_Proud,
                                 gk.estimate_Angry, gk.estimate_Hope,
                                 gk.estimate_Worry, gk.estimate_Excited)

gk.estimate.rap$Estimate = round(gk.estimate.rap$Estimate, 3)
gk.estimate.rap$SE = round(gk.estimate.rap$SE, 3)
gk.estimate.rap$t = round(gk.estimate.rap$t, 3)
gk.estimate.rap$p = round(gk.estimate.rap$p, 3)

write.csv(gk.estimate.rap, file = file.path(output_dir, "11_2_high_interest_potential_mechanisms.csv"), row.names = TRUE)



# 11.2.4 low.interest - potential mechanisms
set.seed(12345)

low.interest$effort = as.numeric(low.interest$effort)
low.interest$trivialization = as.numeric(low.interest$trivialization)
low.interest$`emotion_1` = as.numeric(low.interest$`emotion_1`)
low.interest$`emotion_2` = as.numeric(low.interest$`emotion_2`)
low.interest$`emotion_3` = as.numeric(low.interest$`emotion_3`)
low.interest$`emotion_4` = as.numeric(low.interest$`emotion_4`)
low.interest$`emotion_5` = as.numeric(low.interest$`emotion_5`)
low.interest$`emotion_6` = as.numeric(low.interest$`emotion_6`)

gk.estimate=list()
gk.estimate[[1]]=estimate(rand = low.interest$rand, tr = low.interest$tr, y = low.interest$effort)
gk.estimate[[2]]=estimate(rand = low.interest$rand, tr = low.interest$tr, y = low.interest$trivialization)
gk.estimate[[3]]=estimate(rand = low.interest$rand, tr = low.interest$tr, y = low.interest$`emotion_1`)
gk.estimate[[4]]=estimate(rand = low.interest$rand, tr = low.interest$tr, y = low.interest$`emotion_2`)
gk.estimate[[5]]=estimate(rand = low.interest$rand, tr = low.interest$tr, y = low.interest$`emotion_3`)
gk.estimate[[6]]=estimate(rand = low.interest$rand, tr = low.interest$tr, y = low.interest$`emotion_4`)
gk.estimate[[7]]=estimate(rand = low.interest$rand, tr = low.interest$tr, y = low.interest$`emotion_5`)
gk.estimate[[8]]=estimate(rand = low.interest$rand, tr = low.interest$tr, y = low.interest$`emotion_6`)

list.name=c("Effort","Trivialization",
            "Anxious","Proud",
            "Angry","Hope",
            "Worry","Excited")

names(gk.estimate)=list.name # remember to give every estimate a name

gk.estimate_Effort = as.data.frame(gk.estimate$Effort)
gk.estimate_Trivialization = as.data.frame(gk.estimate$Trivialization)
gk.estimate_Anxious = as.data.frame(gk.estimate$Anxious)
gk.estimate_Proud = as.data.frame(gk.estimate$Proud)
gk.estimate_Angry = as.data.frame(gk.estimate$Angry)
gk.estimate_Hope = as.data.frame(gk.estimate$Hope)
gk.estimate_Worry = as.data.frame(gk.estimate$Worry)
gk.estimate_Excited = as.data.frame(gk.estimate$Excited)

gk.estimate_Effort['name'] = rep("Effort", times=4)
gk.estimate_Trivialization['name'] = rep("Trivialization", times=4)
gk.estimate_Anxious['name'] = rep("Anxious", times=4)
gk.estimate_Proud['name'] = rep("Proud", times=4)
gk.estimate_Angry['name'] = rep("Angry", times=4)
gk.estimate_Hope['name'] = rep("Hope", times=4)
gk.estimate_Worry['name'] = rep("Worry", times=4)
gk.estimate_Excited['name'] = rep("Excited", times=4)

gk.estimate.rap=rbind.data.frame(gk.estimate_Effort, gk.estimate_Trivialization,
                                 gk.estimate_Anxious, gk.estimate_Proud,
                                 gk.estimate_Angry, gk.estimate_Hope,
                                 gk.estimate_Worry, gk.estimate_Excited)

gk.estimate.rap$Estimate = round(gk.estimate.rap$Estimate, 3)
gk.estimate.rap$SE = round(gk.estimate.rap$SE, 3)
gk.estimate.rap$t = round(gk.estimate.rap$t, 3)
gk.estimate.rap$p = round(gk.estimate.rap$p, 3)

write.csv(gk.estimate.rap, file = file.path(output_dir, "11_2_low_interest_potential_mechanisms.csv"), row.names = TRUE)



# 11.3 political news consumption
data$news_consumption = as.numeric(data$news_consumption)
table(data$news_consumption) # 0:973, 0.25:711, 0.5:1628, 0.75:1046, 1:807

data$high_news_consumption=ifelse(std(data$news_consumption)>=0,1,0)
table(data$high_news_consumption) # 0:3312, 1:1853

high.news_consumption <- data[data$high_news_consumption == 1,]
low.news_consumption <- data[data$high_news_consumption == 0,]


# 11.3.1 high.news_consumption
dt_high.news_consumption <- as.data.frame(matrix(nrow=nrow(high.news_consumption)*2,ncol=0)) 
dt_high.news_consumption$central <- append(high.news_consumption$central_trust, high.news_consumption$central_satisfaction)
dt_high.news_consumption$local <- append(high.news_consumption$local_trust, high.news_consumption$local_satisfaction)
dt_high.news_consumption$evaluation <- append(high.news_consumption$evaluation_now, high.news_consumption$evaluation_future)
dt_high.news_consumption$confidence <- append(high.news_consumption$confidence_1, high.news_consumption$confidence_2)
dt_high.news_consumption$efficacy <- append(high.news_consumption$efficacy_1, high.news_consumption$efficacy_2)
dt_high.news_consumption$pride <- append(high.news_consumption$pride_1, high.news_consumption$pride_2)
dt_high.news_consumption$rand <- append(high.news_consumption$rand, high.news_consumption$rand)
dt_high.news_consumption$tr <- append(high.news_consumption$tr, high.news_consumption$tr)

result_dt_high.news_consumption=principal(dt_high.news_consumption[,c(1,2,4,5)], nfactors = 1, residuals = FALSE,rotate="varimax")
result_dt_high.news_consumption
dt_high.news_consumption$pca_score=result_dt_high.news_consumption$score[,1]
dt_high.news_consumption$pca_score

gk_11.3.1.estimate=list()
gk_11.3.1.estimate[[1]] = estimate(rand = dt_high.news_consumption$rand, tr = dt_high.news_consumption$tr, y = dt_high.news_consumption$pca_score)

list_.name=c("PCA Score")
names(gk_11.3.1.estimate)=list_.name
gk_11.3.1.estimate_pca_score = as.data.frame(gk_11.3.1.estimate$`PCA Score`)
gk_11.3.1.estimate_pca_score['name'] = rep("PCA Score", times=4)

gk_11.3.1.estimate_pca_score$Estimate = round(gk_11.3.1.estimate_pca_score$Estimate, 3)
gk_11.3.1.estimate_pca_score$SE = round(gk_11.3.1.estimate_pca_score$SE, 3)
gk_11.3.1.estimate_pca_score$t = round(gk_11.3.1.estimate_pca_score$t, 3)
gk_11.3.1.estimate_pca_score$p = round(gk_11.3.1.estimate_pca_score$p, 3)
write.csv(gk_11.3.1.estimate_pca_score, file = file.path(output_dir, "11_3_pca_on_pro_regime_attitudes_high_news_consumption.csv"), row.names = TRUE)



# 11.3.2 low.news_consumption
dt_low.news_consumption <- as.data.frame(matrix(nrow=nrow(low.news_consumption)*2,ncol=0)) 
dt_low.news_consumption$central <- append(low.news_consumption$central_trust, low.news_consumption$central_satisfaction)
dt_low.news_consumption$local <- append(low.news_consumption$local_trust, low.news_consumption$local_satisfaction)
dt_low.news_consumption$evaluation <- append(low.news_consumption$evaluation_now, low.news_consumption$evaluation_future)
dt_low.news_consumption$confidence <- append(low.news_consumption$confidence_1, low.news_consumption$confidence_2)
dt_low.news_consumption$efficacy <- append(low.news_consumption$efficacy_1, low.news_consumption$efficacy_2)
dt_low.news_consumption$pride <- append(low.news_consumption$pride_1, low.news_consumption$pride_2)
dt_low.news_consumption$rand <- append(low.news_consumption$rand, low.news_consumption$rand)
dt_low.news_consumption$tr <- append(low.news_consumption$tr, low.news_consumption$tr)

result_dt_low.news_consumption=principal(dt_low.news_consumption[,c(1,2,4,5)], nfactors = 1, residuals = FALSE,rotate="varimax")
result_dt_low.news_consumption
dt_low.news_consumption$pca_score=result_dt_low.news_consumption$score[,1]
dt_low.news_consumption$pca_score


gk_11.3.2.estimate=list()
gk_11.3.2.estimate[[1]] = estimate(rand = dt_low.news_consumption$rand, tr = dt_low.news_consumption$tr, y = dt_low.news_consumption$pca_score)

list_.name=c("PCA Score")

names(gk_11.3.2.estimate)=list_.name

gk_11.3.2.estimate_pca_score = as.data.frame(gk_11.3.2.estimate$`PCA Score`)

gk_11.3.2.estimate_pca_score['name'] = rep("PCA Score", times=4)

gk_11.3.2.estimate_pca_score$Estimate = round(gk_11.3.2.estimate_pca_score$Estimate, 3)
gk_11.3.2.estimate_pca_score$SE = round(gk_11.3.2.estimate_pca_score$SE, 3)
gk_11.3.2.estimate_pca_score$t = round(gk_11.3.2.estimate_pca_score$t, 3)
gk_11.3.2.estimate_pca_score$p = round(gk_11.3.2.estimate_pca_score$p, 3)
write.csv(gk_11.3.2.estimate_pca_score, file = file.path(output_dir, "11_3_pca_on_pro_regime_attitudes_low_news_consumption.csv"), row.names = TRUE)



### 12 Figures ###

# 12.1 Propaganda Effects on PCA Scores of Pro-regime Attitudes

header_s1_pca <- read.csv(file.path(data_dir_survey1, "9_pca_on_pro_regime_attitudes.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_s1_pca <- read.csv(file.path(data_dir_survey1, "9_pca_on_pro_regime_attitudes.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_s1_pca) = header_s1_pca
rm(header_s1_pca) 
data_s1_pca <- data_s1_pca[-1, -1]
data_s1_pca$`Surveys (propaganda vs. non-propaganda)` <- 'Survey 1'
data_s1_pca$Estimate <- as.numeric(data_s1_pca$Estimate)
data_s1_pca$SE <- as.numeric(data_s1_pca$SE)
data_s1_pca <- data_s1_pca %>%
  filter(Effect != "naive")

header_s2_pca <- read.csv(file.path(data_dir_survey2, "9_pca_on_pro_regime_attitudes.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_s2_pca <- read.csv(file.path(data_dir_survey2, "9_pca_on_pro_regime_attitudes.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_s2_pca) = header_s2_pca
rm(header_s2_pca) 
data_s2_pca <- data_s2_pca[-1, -1]
data_s2_pca$`Surveys (propaganda vs. non-propaganda)` <- 'Survey 2'
data_s2_pca$Estimate <- as.numeric(data_s2_pca$Estimate)
data_s2_pca$SE <- as.numeric(data_s2_pca$SE)
data_s2_pca <- data_s2_pca %>%
  filter(Effect != "naive")

header_s3_pca <- read.csv(file.path(data_dir_survey3, "9_pca_on_pro_regime_attitudes.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_s3_pca <- read.csv(file.path(data_dir_survey3, "9_pca_on_pro_regime_attitudes.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_s3_pca) = header_s3_pca
rm(header_s3_pca) 
data_s3_pca <- data_s3_pca[-1, -1]
data_s3_pca$`Surveys (propaganda vs. non-propaganda)` <- 'Survey 3'
data_s3_pca$Estimate <- as.numeric(data_s3_pca$Estimate)
data_s3_pca$SE <- as.numeric(data_s3_pca$SE)
data_s3_pca <- data_s3_pca %>%
  filter(Effect != "naive")

header_s4_pca <- read.csv(file.path(data_dir_survey4, "9_pca_on_pro_regime_attitudes.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_s4_pca <- read.csv(file.path(data_dir_survey4, "9_pca_on_pro_regime_attitudes.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_s4_pca) = header_s4_pca
rm(header_s4_pca) 
data_s4_pca <- data_s4_pca[-1, -1]
data_s4_pca$`Surveys (propaganda vs. non-propaganda)` <- 'Survey 4'
data_s4_pca$Estimate <- as.numeric(data_s4_pca$Estimate)
data_s4_pca$SE <- as.numeric(data_s4_pca$SE)
data_s4_pca <- data_s4_pca %>%
  filter(Effect != "naive")

header_s5_pca <- read.csv(file.path(data_dir_survey5, "9_pca_on_pro_regime_attitudes.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_s5_pca <- read.csv(file.path(data_dir_survey5, "9_pca_on_pro_regime_attitudes.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_s5_pca) = header_s5_pca
rm(header_s5_pca) 
data_s5_pca <- data_s5_pca[-1, -1]
data_s5_pca$`Surveys (propaganda vs. non-propaganda)` <- 'Survey 5'
data_s5_pca$Estimate <- as.numeric(data_s5_pca$Estimate)
data_s5_pca$SE <- as.numeric(data_s5_pca$SE)
data_s5_pca <- data_s5_pca %>%
  filter(Effect != "naive")

header_s6_pca <- read.csv(file.path(data_dir_survey6, "9_pca_on_pro_regime_attitudes.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_s6_pca <- read.csv(file.path(data_dir_survey6, "9_pca_on_pro_regime_attitudes.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_s6_pca) = header_s6_pca
rm(header_s6_pca) 
data_s6_pca <- data_s6_pca[-1, -1]
data_s6_pca$`Surveys (propaganda vs. non-propaganda)` <- 'Survey 6'
data_s6_pca$Estimate <- as.numeric(data_s6_pca$Estimate)
data_s6_pca$SE <- as.numeric(data_s6_pca$SE)
data_s6_pca <- data_s6_pca %>%
  filter(Effect != "naive")

header_combine_pca <- read.csv(file.path(output_dir, "9_pca_on_pro_regime_attitudes.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_combine_pca <- read.csv(file.path(output_dir, "9_pca_on_pro_regime_attitudes.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_combine_pca) = header_combine_pca
rm(header_combine_pca) 
data_combine_pca <- data_combine_pca[-1, -1]
data_combine_pca$`Surveys (propaganda vs. non-propaganda)` <- 'Combined'
data_combine_pca$Estimate <- as.numeric(data_combine_pca$Estimate)
data_combine_pca$SE <- as.numeric(data_combine_pca$SE)
data_combine_pca <- data_combine_pca %>%
  filter(Effect != "naive")

data_pca = rbind.data.frame(data_combine_pca, data_s1_pca, data_s2_pca, data_s3_pca, data_s4_pca, data_s5_pca, data_s6_pca)

data_pca <- data_pca %>%
  mutate(Effect = ifelse(Effect == "t", "ATE", Effect)) %>%
  mutate(Effect = ifelse(Effect == "t_s", "Selector", Effect)) %>%
  mutate(Effect = ifelse(Effect == "t_n", "Non-selector", Effect))

pd <- position_dodge(0.5)

figure_pca <- ggplot(data_pca, aes(x=Effect, y=Estimate, group=`Surveys (propaganda vs. non-propaganda)`)) + 
  geom_errorbar(aes(ymin=Estimate-1.96*SE, ymax=Estimate+1.96*SE), width=0, position=pd) +
  geom_point(aes(shape=`Surveys (propaganda vs. non-propaganda)`), position=pd)+scale_shape_manual(values=seq(0,7))+
  xlab("Group") +
  ylab("PCA Scores of Pro-regime Attitudes") +
  expand_limits(y=0) +                        # Expand y range
  ylim(-1.2,2.4) +         # Set tick every 4
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_hline(yintercept=0,color="red",linetype="dotted")

ggsave(filename = file.path(output_dir, "12_1_Propaganda Effects on PCA Scores of Pro-regime Attitudes.png"), plot = figure_pca, width = 9, height = 3, dpi = 300)



# 12.2 Marginal Effects of Independent Variables on Reading Preferences (Combined: propaganda vs. non-propaganda) 

header_combine_marginal_effects <- read.csv(file.path(output_dir, "8_marginal_effects.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_combine_marginal_effects <- read.csv(file.path(output_dir, "8_marginal_effects.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_combine_marginal_effects) = header_combine_marginal_effects
rm(header_combine_marginal_effects) 
data_combine_marginal_effects <- data_combine_marginal_effects[-c(1, 2), ]

data_combine_marginal_effects$Estimate = as.numeric(data_combine_marginal_effects$Estimate)
data_combine_marginal_effects$`Std. Error` = as.numeric(data_combine_marginal_effects$`Std. Error`)

data_combine_marginal_effects <- data_combine_marginal_effects %>%
  mutate(`NA` = ifelse(`NA` == "age", "Age", `NA`)) %>%
  mutate(`NA` = ifelse(`NA` == "education", "Education", `NA`)) %>%
  mutate(`NA` = ifelse(`NA` == "ethnicity", "Ethnicity", `NA`)) %>%
  mutate(`NA` = ifelse(`NA` == "gender", "Gender", `NA`)) %>%
  mutate(`NA` = ifelse(`NA` == "income", "Income", `NA`)) %>%
  mutate(`NA` = ifelse(`NA` == "interest", "Interest", `NA`)) %>%
  mutate(`NA` = ifelse(`NA` == "knowledge", "Knowledge", `NA`)) %>%
  mutate(`NA` = ifelse(`NA` == "major", "Major", `NA`)) %>%
  mutate(`NA` = ifelse(`NA` == "marriage", "Marriage", `NA`)) %>%
  mutate(`NA` = ifelse(`NA` == "news_consumption", "News consumption", `NA`)) %>%
  mutate(`NA` = ifelse(`NA` == "occupation", "Occupation", `NA`)) %>%
  mutate(`NA` = ifelse(`NA` == "party", "Party", `NA`)) %>%
  mutate(`NA` = ifelse(`NA` == "urban", "Urban", `NA`))

figure_marginal_effects <- ggplot(data_combine_marginal_effects, aes(x=`NA`, y=Estimate)) + 
  geom_errorbar(aes(ymin=Estimate-1.96*`Std. Error`, ymax=Estimate+1.96*`Std. Error`), width=0, position=pd) +
  geom_point(position=pd) + scale_shape_manual(values=seq(0,7))+
  xlab("") +
  ylab("Mean") +
  expand_limits(y=0) +                        # Expand y range
  ylim(-0.6,1.1) +         # Set tick every 4
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_hline(yintercept=0,color="red",linetype="dotted")

ggsave(filename = file.path(output_dir, "12_2_Marginal Effects of Independent Variables on Reading Preferences (Combined: propaganda vs. non-propaganda).png"), plot = figure_marginal_effects, width = 12, height = 3, dpi = 300)



# 12.3 Potential Mechanisms (Effort justification)

header_s1_effort <- read.csv(file.path(data_dir_survey1, "10_potential_mechanisms.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_s1_effort <- read.csv(file.path(data_dir_survey1, "10_potential_mechanisms.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_s1_effort) = header_s1_effort
rm(header_s1_effort) 
data_s1_effort <- data_s1_effort[-1, -1]
data_s1_effort <- data_s1_effort %>%
  filter(name == "Effort") %>%
  filter(Effect != "naive")
data_s1_effort$`Surveys (propaganda vs. non-propaganda)` <- 'Survey 1'
data_s1_effort$Estimate <- as.numeric(data_s1_effort$Estimate)
data_s1_effort$SE <- as.numeric(data_s1_effort$SE)

header_s2_effort <- read.csv(file.path(data_dir_survey2, "10_potential_mechanisms.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_s2_effort <- read.csv(file.path(data_dir_survey2, "10_potential_mechanisms.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_s2_effort) = header_s2_effort
rm(header_s2_effort) 
data_s2_effort <- data_s2_effort[-1, -1]
data_s2_effort <- data_s2_effort %>%
  filter(name == "Effort") %>%
  filter(Effect != "naive")
data_s2_effort$`Surveys (propaganda vs. non-propaganda)` <- 'Survey 2'
data_s2_effort$Estimate <- as.numeric(data_s2_effort$Estimate)
data_s2_effort$SE <- as.numeric(data_s2_effort$SE)

header_s3_effort <- read.csv(file.path(data_dir_survey3, "10_potential_mechanisms.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_s3_effort <- read.csv(file.path(data_dir_survey3, "10_potential_mechanisms.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_s3_effort) = header_s3_effort
rm(header_s3_effort) 
data_s3_effort <- data_s3_effort[-1, -1]
data_s3_effort <- data_s3_effort %>%
  filter(name == "Effort") %>%
  filter(Effect != "naive")
data_s3_effort$`Surveys (propaganda vs. non-propaganda)` <- 'Survey 3'
data_s3_effort$Estimate <- as.numeric(data_s3_effort$Estimate)
data_s3_effort$SE <- as.numeric(data_s3_effort$SE)

header_s4_effort <- read.csv(file.path(data_dir_survey4, "10_potential_mechanisms.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_s4_effort <- read.csv(file.path(data_dir_survey4, "10_potential_mechanisms.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_s4_effort) = header_s4_effort
rm(header_s4_effort) 
data_s4_effort <- data_s4_effort[-1, -1]
data_s4_effort <- data_s4_effort %>%
  filter(name == "Effort") %>%
  filter(Effect != "naive")
data_s4_effort$`Surveys (propaganda vs. non-propaganda)` <- 'Survey 4'
data_s4_effort$Estimate <- as.numeric(data_s4_effort$Estimate)
data_s4_effort$SE <- as.numeric(data_s4_effort$SE)

header_s5_effort <- read.csv(file.path(data_dir_survey5, "10_potential_mechanisms.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_s5_effort <- read.csv(file.path(data_dir_survey5, "10_potential_mechanisms.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_s5_effort) = header_s5_effort
rm(header_s5_effort) 
data_s5_effort <- data_s5_effort[-1, -1]
data_s5_effort <- data_s5_effort %>%
  filter(name == "Effort") %>%
  filter(Effect != "naive")
data_s5_effort$`Surveys (propaganda vs. non-propaganda)` <- 'Survey 5'
data_s5_effort$Estimate <- as.numeric(data_s5_effort$Estimate)
data_s5_effort$SE <- as.numeric(data_s5_effort$SE)

header_s6_effort <- read.csv(file.path(data_dir_survey6, "10_potential_mechanisms.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_s6_effort <- read.csv(file.path(data_dir_survey6, "10_potential_mechanisms.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_s6_effort) = header_s6_effort
rm(header_s6_effort) 
data_s6_effort <- data_s6_effort[-1, -1]
data_s6_effort <- data_s6_effort %>%
  filter(name == "Effort") %>%
  filter(Effect != "naive")
data_s6_effort$`Surveys (propaganda vs. non-propaganda)` <- 'Survey 6'
data_s6_effort$Estimate <- as.numeric(data_s6_effort$Estimate)
data_s6_effort$SE <- as.numeric(data_s6_effort$SE)

header_combine_effort <- read.csv(file.path(data_dir_survey6, "10_potential_mechanisms.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_combine_effort <- read.csv(file.path(data_dir_survey6, "10_potential_mechanisms.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_combine_effort) = header_combine_effort
rm(header_combine_effort) 
data_combine_effort <- data_combine_effort[-1, -1]
data_combine_effort <- data_combine_effort %>%
  filter(name == "Effort") %>%
  filter(Effect != "naive")
data_combine_effort$`Surveys (propaganda vs. non-propaganda)` <- 'Combined'
data_combine_effort$Estimate <- as.numeric(data_combine_effort$Estimate)
data_combine_effort$SE <- as.numeric(data_combine_effort$SE)

data_effort = rbind.data.frame(data_combine_effort, data_s1_effort, data_s2_effort, data_s3_effort, data_s4_effort, data_s5_effort, data_s6_effort)

data_effort <- data_effort %>%
  mutate(Effect = ifelse(Effect == "t", "ATE", Effect)) %>%
  mutate(Effect = ifelse(Effect == "t_s", "Selector", Effect)) %>%
  mutate(Effect = ifelse(Effect == "t_n", "Non-selector", Effect))

pd <- position_dodge(0.5)

figure_effort <- ggplot(data_effort, aes(x=Effect, y=Estimate, group=`Surveys (propaganda vs. non-propaganda)`)) + 
  geom_errorbar(aes(ymin=Estimate-1.96*SE, ymax=Estimate+1.96*SE), width=0, position=pd) +
  geom_point(aes(shape=`Surveys (propaganda vs. non-propaganda)`), position=pd)+scale_shape_manual(values=seq(0,7))+
  xlab("Group") +
  ylab("Mean") +
  expand_limits(y=0) +                        # Expand y range
  ylim(-1.2,2.4) +         # Set tick every 4
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_hline(yintercept=0,color="red",linetype="dotted")

ggsave(filename = file.path(output_dir, "12_3_Potential Mechanisms (Effort justification).png"), plot = figure_effort, width = 9, height = 3, dpi = 300)



# 12.4 Potential Mechanisms (Trivialization)

header_s1_trivialization <- read.csv(file.path(data_dir_survey1, "10_potential_mechanisms.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_s1_trivialization <- read.csv(file.path(data_dir_survey1, "10_potential_mechanisms.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_s1_trivialization) = header_s1_trivialization
rm(header_s1_trivialization) 
data_s1_trivialization <- data_s1_trivialization[-1, -1]
data_s1_trivialization <- data_s1_trivialization %>%
  filter(name == "Trivialization") %>%
  filter(Effect != "naive")
data_s1_trivialization$`Surveys (propaganda vs. non-propaganda)` <- 'Survey 1'
data_s1_trivialization$Estimate <- as.numeric(data_s1_trivialization$Estimate)
data_s1_trivialization$SE <- as.numeric(data_s1_trivialization$SE)

header_s2_trivialization <- read.csv(file.path(data_dir_survey2, "10_potential_mechanisms.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_s2_trivialization <- read.csv(file.path(data_dir_survey2, "10_potential_mechanisms.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_s2_trivialization) = header_s2_trivialization
rm(header_s2_trivialization) 
data_s2_trivialization <- data_s2_trivialization[-1, -1]
data_s2_trivialization <- data_s2_trivialization %>%
  filter(name == "Trivialization") %>%
  filter(Effect != "naive")
data_s2_trivialization$`Surveys (propaganda vs. non-propaganda)` <- 'Survey 2'
data_s2_trivialization$Estimate <- as.numeric(data_s2_trivialization$Estimate)
data_s2_trivialization$SE <- as.numeric(data_s2_trivialization$SE)

header_s3_trivialization <- read.csv(file.path(data_dir_survey3, "10_potential_mechanisms.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_s3_trivialization <- read.csv(file.path(data_dir_survey3, "10_potential_mechanisms.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_s3_trivialization) = header_s3_trivialization
rm(header_s3_trivialization) 
data_s3_trivialization <- data_s3_trivialization[-1, -1]
data_s3_trivialization <- data_s3_trivialization %>%
  filter(name == "Trivialization") %>%
  filter(Effect != "naive")
data_s3_trivialization$`Surveys (propaganda vs. non-propaganda)` <- 'Survey 3'
data_s3_trivialization$Estimate <- as.numeric(data_s3_trivialization$Estimate)
data_s3_trivialization$SE <- as.numeric(data_s3_trivialization$SE)

header_s4_trivialization <- read.csv(file.path(data_dir_survey4, "10_potential_mechanisms.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_s4_trivialization <- read.csv(file.path(data_dir_survey4, "10_potential_mechanisms.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_s4_trivialization) = header_s4_trivialization
rm(header_s4_trivialization) 
data_s4_trivialization <- data_s4_trivialization[-1, -1]
data_s4_trivialization <- data_s4_trivialization %>%
  filter(name == "Trivialization") %>%
  filter(Effect != "naive")
data_s4_trivialization$`Surveys (propaganda vs. non-propaganda)` <- 'Survey 4'
data_s4_trivialization$Estimate <- as.numeric(data_s4_trivialization$Estimate)
data_s4_trivialization$SE <- as.numeric(data_s4_trivialization$SE)

header_s5_trivialization <- read.csv(file.path(data_dir_survey5, "10_potential_mechanisms.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_s5_trivialization <- read.csv(file.path(data_dir_survey5, "10_potential_mechanisms.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_s5_trivialization) = header_s5_trivialization
rm(header_s5_trivialization) 
data_s5_trivialization <- data_s5_trivialization[-1, -1]
data_s5_trivialization <- data_s5_trivialization %>%
  filter(name == "Trivialization") %>%
  filter(Effect != "naive")
data_s5_trivialization$`Surveys (propaganda vs. non-propaganda)` <- 'Survey 5'
data_s5_trivialization$Estimate <- as.numeric(data_s5_trivialization$Estimate)
data_s5_trivialization$SE <- as.numeric(data_s5_trivialization$SE)

header_s6_trivialization <- read.csv(file.path(data_dir_survey6, "10_potential_mechanisms.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_s6_trivialization <- read.csv(file.path(data_dir_survey6, "10_potential_mechanisms.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_s6_trivialization) = header_s6_trivialization
rm(header_s6_trivialization) 
data_s6_trivialization <- data_s6_trivialization[-1, -1]
data_s6_trivialization <- data_s6_trivialization %>%
  filter(name == "Trivialization") %>%
  filter(Effect != "naive")
data_s6_trivialization$`Surveys (propaganda vs. non-propaganda)` <- 'Survey 6'
data_s6_trivialization$Estimate <- as.numeric(data_s6_trivialization$Estimate)
data_s6_trivialization$SE <- as.numeric(data_s6_trivialization$SE)

header_combine_trivialization <- read.csv(file.path(data_dir_survey6, "10_potential_mechanisms.csv"), header = FALSE, nrows = 1, as.is = TRUE)
data_combine_trivialization <- read.csv(file.path(data_dir_survey6, "10_potential_mechanisms.csv"), header= FALSE , na=c("", "NA")) 
colnames(data_combine_trivialization) = header_combine_trivialization
rm(header_combine_trivialization) 
data_combine_trivialization <- data_combine_trivialization[-1, -1]
data_combine_trivialization <- data_combine_trivialization %>%
  filter(name == "Trivialization") %>%
  filter(Effect != "naive")
data_combine_trivialization$`Surveys (propaganda vs. non-propaganda)` <- 'Combined'
data_combine_trivialization$Estimate <- as.numeric(data_combine_trivialization$Estimate)
data_combine_trivialization$SE <- as.numeric(data_combine_trivialization$SE)

data_trivialization = rbind.data.frame(data_combine_trivialization, data_s1_trivialization, data_s2_trivialization, data_s3_trivialization, data_s4_trivialization, data_s5_trivialization, data_s6_trivialization)

data_trivialization <- data_trivialization %>%
  mutate(Effect = ifelse(Effect == "t", "ATE", Effect)) %>%
  mutate(Effect = ifelse(Effect == "t_s", "Selector", Effect)) %>%
  mutate(Effect = ifelse(Effect == "t_n", "Non-selector", Effect))

pd <- position_dodge(0.5)

figure_trivialization <- ggplot(data_trivialization, aes(x=Effect, y=Estimate, group=`Surveys (propaganda vs. non-propaganda)`)) + 
  geom_errorbar(aes(ymin=Estimate-1.96*SE, ymax=Estimate+1.96*SE), width=0, position=pd) +
  geom_point(aes(shape=`Surveys (propaganda vs. non-propaganda)`), position=pd)+scale_shape_manual(values=seq(0,7))+
  xlab("Group") +
  ylab("Mean") +
  expand_limits(y=0) +                        # Expand y range
  ylim(-1.2,2.4) +         # Set tick every 4
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_hline(yintercept=0,color="red",linetype="dotted")

ggsave(filename = file.path(output_dir, "12_4_Potential Mechanisms (Trivialization).png"), plot = figure_trivialization, width = 9, height = 3, dpi = 300)


